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Characterisation of global flow and local fluctuations in 3D 
SPH simulations of protoplanetary discs 
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ABSTRACT 

A complete and detailed knowledge of the structure of the gaseous component in pro- 
toplanetary discs is essential to the study of dust evolution during the early phases of 
pre-planetesimal formation. The aim of this paper is to determine if three-dimensional 
accretion discs simulated by the Smoothed Particle Hydrodynamics (SPH) method can 
reproduce the observational data now available and the expected turbulent nature of 
protoplanetary discs. The investigation is carried out by setting up a suite of diagnostic 
tools specifically designed to characterise both the global flow and the fluctuations of 
the gaseous disc. The main result concerns the role of the artificial viscosity implemen- 
tation in the SPH method: in addition to the already known ability of SPH artificial 
viscosity to mimic a physical-like viscosity under specific conditions, we show how the 
same artificial viscosity prescription behaves like an implicit turbulence model. In fact, 
we identify a threshold for the parameters in the standard artificial viscosity above 
which SPH disc models present a cascade in the power spectrum of velocity fluctua- 
tions, turbulent diffusion and a mass accretion rate of the same order of magnitude as 
measured in observations. Furthermore, the turbulence properties observed locally in 
SPH disc models are accompanied by meridional circulation in the global flow of the 
gas, proving that the two mechanisms can coexist. 

Key words: accretion, accretion discs - hydrodynamics - turbulence - methods: 
numerical - protoplanetary discs. 



1 INTRODUCTION 

Protoplanetary discs (abbreviated PPD) are discs composed 
mainly of gas and dust in quasi-Keplerian rotation around 
young stars. They are believed to be the birth place of plan- 
ets. The theory of planet formation is complex because sev- 
eral scales (from /im to hundreds of thousands of km) , sev- 
eral forces (e.g. electrostatic, magnetic, gravitational, drag) 
and several processes (e.g. turbulence, chemical and thermo- 
dynamical transformations, instabilities) are involved. One 
of the more poorly understood stages is the evolution of /im- 
size dust grains into km-size objects, called planetesimals. 

Among the several mechanisms the dust is subject to 
(e.g. radial drift and vertical settling), possible large scale 
motion and turbulence are expected to be of relevant im- 
portance. Large scale motion can lead dust to travel to very 
different locations in the disc. Turbulence can have two com- 
peting effects mediated by gas drag: stirring up and diffusing 
dust particles, impeding their agglomeration, or trapping 
them inside eddies, favoring their agglomeration (see e.i 
ICuzzi et al.lll993l : ICarballido eralll2008l : ICuzzi et allboOS 

Here we focus on the global flow and on the expected 



turbulent behaviour of only the gaseous component of PPD, 
which represents the medium in which dust evolves. 

In the literature, two large classes of models of turbu- 
lent discs are widely used: continuous and discrete models. 
Continuous models of viscous discs are directly derived from 
the Navier-Stokes equations in presence of the gravitational 
potential of the central star, therefore they include the in- 
gredients of physical viscosity and star gravity. Among the 
most pop ular models of this type we find the analytic ID 
models bv lPringld jlQSlI ') and lLynden-Bell fc Pringlel(|l974 ). 
Since the physical molecular visco sity of the gas in PPD is 
very low (see e.g. I Armitagd 120071 ). high Reynolds numbers 
are expected and therefore the gas is believed to be turbu- 
lent. The most commonly used model for such discs is the 
IShakura fc Sunvaevl (jl973l l model (hereafter SS73), which 
is a Prandtl model for turbulence: the Navier-Stokes equa- 
tions are averaged and the Reynolds stress tensor is mod- 
eled by a viscosity called turbulent viscosity. The resulting 
total viscosity is the sum of the molecular and the turbulent 
viscosities, the former being negligible in PPD. Continuous 
models are usually applied to turbulent discs whose source 
of turbulence is not known or not explicitly declared; in such 



S. E. Arena and J.-F. Gonzalez 



cases the|SS73 relation vt — aaaCaH is used to express the 
turbulent viscosity parameter. The dimensionless parameter 
Oss, usually taken as a constant, collects all the ignorance 
of the source of turbulence, Cg is the gas sound speed and 
H the scale height of the disc. Generally, ass-discs are good 
models for turbulence and capture its main effects. An ex- 
ample is given by the t wo-dimensional models calculated by 
iTakeuchi fc LinI (|2002l . hereafter TL02) that show how the 
ass prescription is able to describe the mechanism of merid- 
ional circulation in discs. However, turbulent fluctuations 
are not directly reproduced in such disc models, only their 
effect on the dynamics of the global discs can be studied. In 
addition, if necessary, turbulent diffusion has to be added to 
the basic equations. 

A different line of research starts from a proposed ex- 
plicit hypothesis concerning the source of turbulence and 
derive the corresponding disc evolution. Given the com- 
plex dynamics, numerical simulations are often required, 
therefore these consist of mainly discrete models. Turbu- 
lent fluctuations can be resolved down to the resolution 
scale, which limits the size of the smallest structures. A first 
example of such models are discs produced by Magneto- 
hydrodynamics (MHD) simulations where the source of tur- 
bulence is the Magneto- Rotational Instability (MRI) (e.g. 
iFromang fc NelsonlbOOa ). The equations at the base of such 
models are the Euler equations in presence of an adequate 
magnetic field and of the gravitational potential of the cen- 
tral star. Therefore they include the ingredients of magnetic 
field and star gravity, however physical viscosity (included 
in continuous models) is not present. In addition, turbulent 
diffusion (not added a priori) has been found in the out- 
come of simulations. A second example of discrete models 
are systems produced by simulations of self-gravitating discs 
where the so urce of turbulenc e is the Gravitational Instabil- 
ity (GI) fe.g lRice et al.ll2005l ). The equations at the base of 
such models are the Euler equations in presence of both self- 
gravity and the gravitational potential of the central star. 
In addition, the numerical scheme is completed by an arti- 
ficial viscosity term. Therefore they include the ingredients 
of self-gravity, star gravity and a 'physical-like' viscosity. 

The nature of the source of turbulence in accretion discs 
and particularly in PPD is still an open issue. Even in the ab- 
sence of magnetic fields and self-gravity, other mechanisms 
can be sources of turbulence, like interactions with external 
objects or pure hydrodyna mic and thermal inst abilities, such 
as the Rossby ins tabilitv (ILovelace et al. 1199911 or the baro- 



clinic instability (jKlahr fc Bodenheimei 



20031 ). The possi- 



bility to have a hydrodynamic instability in accretion discs 
is s till debated. In fact, due to the R aylcigh criterion (see 
e.g. lBalbus fc Hawlevlll99ll : lArmitage|[2007; 'l. Keplerian discs 
are believed to be stable with respect to hydrodynamic in- 
stabilities. However this is proved only for ID discs, if the 
third dimension is considered the behaviour could change 
and there may be the possibility to trigger instabilities (see 
e.g. iNelson et ai]|2012l i. 

As done in most continuous analytic models, we pre- 
fer to model the effects of turbulence (and not its cause) 
by the use of an effective viscosity. Here we consider three- 
dimensional accretion discs modeled by means of standard 
Smoothed Parti cle Hydrodynamics (SPH, for a review see 
iMonaghanI |2005| ) . which we call SPH disc models, they are 
discrete models, such as those of the second class. The SPH 



scheme includes an artificial viscosity term (a Von Neumann- 
like viscosity) that has been introduce d to treat shocks cor- 
rectly (see e.g. lCuUen fc Dehnenll20ld V Far from the shock, 
the artificial viscosity is usually turned off by a switch. Dif- 
ferent fiavours of the artificial viscosity term are discussed 
in Sect. O 

The aim of the present paper is to answer the following 
question: can SPH disc models correctly reproduce both the 
observed properties of PPD and the expected effect of tur- 
bulence? In order to clarify these two points we present a 
detailed characterisation of the properties of the global flow 
and particularly of the density and velocity fluctuations of 
SPH disc models. 

The global fiow of SPH models gives us information 
about the average properties of the gaseous disc, while the 
divergence from the average values defines the SPH fluctu- 
ations associated to the considered field (e.g. velocity and 
density fields). In particular, SPH fluctuations result from 
the combination of standard numerical noise and physical 
fluctuations. The former, present in all numerical schemes, 
is due in the case of the SPH method both to the discreti- 
sation of the continuum equations and to the SPH approx- 
imations; it depends on the number of particles and on the 
SPH kernel used in the simulations (see Sect. 12.1]) . The lat- 
ter is present if the gas flow is turbulent; it depends on the 
Reynolds number of the flow which in turn is related to the 
physical viscosity of the gas. Both the numerical noise and 
physical fluctuations are related to the artiflcial viscosity 
term in the specific way we will show in this paper. 

Turbulence modeling with the SP H method has rec ently 
become an active field of research (see lMonagharJl2005l . and 
references therein). The first effort of im plementing tur- 
bulen ce models in the SPH equations (e.g. IVioleau fc Issal 
|2007|) is now moving in the direction of quantifying the 
ability of the SPH method to i ntrinsically repro duce tur- 
bulence (e.g. lEUero et all l2010l : [Monagha3[20lJ). A com- 
parison between SPH a nd grid methods is presented in 
IPrice fc FederrathI (|2010|) and a detailed and clear relation 
between tur bulence and resolution in the SPH method is 
explained in IPricd (|2012r ). In those studies, isotropic ho- 
mogeneous turbulence in 2D boxes with periodic boundary 
conditions and with negligible gravity effects are considered. 
However, the case of anisotropic and inhomogeneous turbu- 
lence with free boundary condition such as in 3D gaseous 
discs, where also the action of gravity from the cen tral star 
is rel e vant, has been address ed only in a few works (|Murravl 
119961 : iLodato fc Pried [2010|) , but without considering the 
global flow and the statistics of the velocity and density field 
and the match with the available observations for PPD. The 
present work aims to contribute in this direction. 

In Sect. [2] we describe the main features of the SPH 
code we have used to model the accretion discs considered in 
the present work, the adopted reference disc model and the 
simulations of the different models we performed. In Sect. [3] 
we present the diagnostics applied to characterise both the 
global gas flow and the fluctuations of the velocity fleld of the 
gas, along with the known reference values for PPD. We then 
show the results of the applied diagnostics respectively for 
the global flow (Sect.|4]), the magnitude of SPH fluctuations 
(Sect. O and the structure of SPH fluctuations (Sect. [S]). 
The results are discussed in Sect. [7] where we also draw our 
conclusions. 
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2 SPH MODELS OF ACCRETION DISCS 

In this work we focus on a PPD described by given analytic 
initial density and locally isothermal sound speed profiles. 
The system is not homogeneous nor isotropic. We calculate 
the quasi-stationary state of such a disc using different val- 
ues of the three numerical parameters A'' (number of parti- 
cles) and a and /3 (artificial viscosity parameters) controlling 
the SPH fluctuations, in order to characterise their effects. 
In this section we present the code and the adopted 
units, the reference disc model and the sets of simulations 
we performed. 



length. The gas is described by a locally isothermal equa- 
tion of state. The artificial viscosity term Ylab is described 
in Sect. [2Xn 



2.1.1 The artificial viscosity term 

Two different implementations of the artificial viscos- 
ity term are considered. The f i rst on e, originally intro- 
duced by iMonaghan fc Gingoldl (|l983h (hereafter MG83 ) 
and subsequently refined by Lattanzio et al.l ([1985|) ; 
iMonaghan fc Lattanziol (|l985l ) ; iMonaghanI l| 19921 ') , is defined 
by: 



2.1 The SPH code 

We use the two-p h ase SPH code described in 
iBarriere-Fouchet et al.l ((20051). The two phases repre- 
sent gas and dust that interact via aerodynamic drag. 
The gas is described by the Euler equations and artificial 
viscosity. 

Here we consider only the gas phase because we are 
interested in the characterisation of the global gas flow and 
its velocity and density fluctuations. Thus, the momentum 
conservation equation takes the form: 



GA'L 
~{r. + ehf'"^' 

where the usual SPH approximation of replacing integrals 
with sums over a finite number of particles A'^ has been per- 
formed. The term Va is the velocity of SPH particle a, t 
the time, m the SPH particle mass, P the pressure, p the 
mass density, Va the position of particle a with respect to 
the central star of mass M, , e = 0.1 a parameter used to 
prevent singularities and W{r,h) the SPH smoothing ker- 
nel, with r = Itq — rfc] the distance between the particle a 
and its neighbour h, and h the smoothing length. Here we 
use a cubic-spline kernel truncated at 2h: 



( 1 



W{q) 



hd 



{2-qf~{l~qf 0<g<l; 







1 < g < 2; ' (2) 



where q — r/h, d is the dimension of the simulation (in this 
work d = 3) and a = [2/3, 10/77r, I/tt] is the normalization 
constant respectively in 1, 2 and 3 dimensions. The effect of 
different kernels will be addressed in a future work. 

The smoothing length h is variable and is derived from 
the density p: 



h 



V 



P 



1/3 



(3) 



with rj = 1.14. This choice of 77 guarantees a roughly 
constant number of neighbours (for 3D simulations: 
A^ncigh ~ 50) . The smoothing length defines the resolution of 
the SPH simulations: for a larger number of particles N the 
target number of neighbours is reached inside a smaller vol- 
ume (of radius 2h from Eq.[2]), implying a smaller smoothing 



n. 



P^J 





-acijfi.j + Pfiij) if 
if 



rij < 
fij > 



(4) 



where rij and Vij are respectively the relative distance and 
relative velocity of particles i and j. The overlined quantities 
are averages between particle i and its neighbouring particle 

j: Cij = {a + Cj)/2, -pij = {pi + pj)/2, hij = {hi + hj)/2. 
Finally, 



Pij — 



IliiVi 



(5) 



with e — 10^ . Different combinations of the two artificial 
viscosity parameters have been used so far in different ap- 
plications, in particular the combin ation (a,p) = (1 , 2) ha s 
been claimed to give good results (|Monaghanlll992l . |2005| ). 
We will refer to a = 1 and f3 — 2 as 'standard' values for 
the artificial viscosity parameters. 

The second implementation of artificial viscosity (here - 
afte r identified as Mu96) is that adopted bv iMurravl (|l996r i 
and iLodato fc Pried (|201G| ') : 



CUCij pij 



(6) 



In contrast to the lMG83l implementation, here the artificial 
viscosi ty is applied to all particles (the switch present in the 
IMG83| version is not applied) and the /3 t erm is set to zero. 
Therefore, discs with the standard IMG83| artificial viscosity 
implementation are naturally half as viscous as those with 
the lMu96l implementation. 



2.1.2 The link to physical viscosity 

The reason of the IMu96| choice is to better match the con- 
ditions under which the artificial viscos i ty can represent a 



physical viscosity. In fact. lMeglicki et al.l ([1993) showed that 
in the 3D continuum limit the IMu96| artificial viscosity has 
the form of a physical (shear and bulk) viscous force. The 
shear term is equivalent to an Oss viscosity (that we call 
Ocont, where the subscript 'cont' refers to continuum) given 
by: 



,„, 1 {h)e.iR) 

«cont,Mu96(.-ftJ ~ -T^Ct- 



10 H{R) 



(7) 



where the average on h is taken along both the azimuthal 
and vertical directions and the semi-thickness H of the disc 
is defined in Sect. 12.21 (note that this expression is valid 
for the cubic spline kernel). However, Qcont depends on the 



S. E. Arena and J.-F. Gonzalez 



reso lution, because of the presence of the smoothing length 
fsee lLodato fc PricdlJOlol ). 

The SPH equations with the iMuQa artificial viscosity 
and with a very large number of particles (high resolution) 
are therefore equivalent to those of a viscous fluid with an 
Qaa given by Eq. [7] Thus, artificial viscosity can be used to 
control the effective viscosity in the disc and it is expected 
to be responsible of physical fluctuations of the flow when 
the gas is in the turbulent regime. 

The 3D continuum limit of the lMGSa artificial vi scosity 
is much more complex than the one presented for the lMuQq 
case and at the moment a simple rel ation such as that in 
Eq.[7]is only available for the a term (JMeru fc Batell2012l ): 



,„, 1 {h)B.{R) 

acont,MG83(-KJ = —a jj,^s ■ 



(8) 



The factor of two between the numerical coefficients of Eqs.[7] 
and [8] naturally arises from the fact that the lMGSa artificial 
viscosity is applied only to half the particles becaus e of the 
presence of the switch (see Eq. [4l l. Since the IMuOb formula 
(developed much earlier than the IMG83| formula) is usually 
adopted to estimate the effective viscosity in SPH discs inde- 
pendently of the particular AV implementation, SPH discs 
are generally less viscous than what has been considered 
so far. In Sect. 14.41 we present a method to determine the 
effective viscosity of an SPH disc with a generic artific ial vis- 
cosity implementation. It is tested using Eq. [7]in the IMu96| 
case. It can therefore be used for sam pling nu merically the 
Qcont.MGsa-ct relation in the full (/3 7^ 0)|MG8^ case, through 
high-resolution simulations. 

We note that in the majority of SPH simulations the use 
of the artificial viscosity is still preferred to the direct imple- 
mentation of the viscous stress tensor to model a pure shear 
viscosity. The main reason is that the latter meth od does 
not conserve the total angu lar momentum exactly (jSchafen 
I2OO5I : iLodato fc Pried I2OICI ) as the artificial viscosity term 
does. Secondarily, it is more computationally demanding be- 
cause of the presence of the second spatial derivative of the 
velocities. 



is the semi-thickness of the disc, related to the sound speed 
Cs and the angular velocity Q. The disc is locally isothermal 
with a temperature radial profile given by 



with g > 0, which leads to the sound speed profile 



Cs{R) = CsO 



-q/2 



(11) 



(12) 



Note that the sound speed coefficient Cso and the sound 
speed exponent q/2 determine respectively the semi- 
thickness of the disc and its radial dependence: 



H(i^)^i/o(| 



(13) 



with Ho = CsqRq /VGM-i,. The resulting radial profile of 
the surface density is 



E(i?) = Eo 



(14) 



with E() = \/^poHo the surface density at Ro and 
p = s+(g-3)/2. 

When only the gas phase is considered, all models are 
self-similar and different physical scales correspond to the 
same dimensionless model. Therefore, in the following, re- 
sults are mainly expressed in code units. 

The reference values we adopt in the following are 
(P) 9) = (3/2, 3/4) and _Ro = 100 au. At this location, the 
disc is slightly flared with Hq/Ro = 0.05, Eo « 4.58 kgm"^ 
and CsO ~ 149 ms~^ {To « 6 K). For reference, the cor- 
responding value at 1 au are: Eq « 4580 kgm~^ and 
To ^ 198 K. 

We follow the evolution of the reference disc model, 
after pressure equilibrium has been reached, up to time 
t = 100 in code units (corresponding to 15.9 orbits at 100 
au). All figures are plotted for that time, unless stated oth- 



2.1.3 The units 

The internal units of the code are chosen to be 1 M© for 
mass, 100 au (Astronomical Units) for length and to give a 
gravitational constant G = 1. With these values the time 
unit is therefore 10"^/27r yr. 



2.2 The reference disc model 

We consider a typical T Tauri disc of mass A/disc ~ 0.01 M* 
orbiting around a one solar mass star (M* = 1 M0). It 
extends from Rin = 20 to -Rout ~ 400 au and is characterised 
by a density profile given, in cylindrical coordinates, by 



p{R,z) 



Po 



exp 



2m ' 



(9) 



(expansion at small z/H, see lLaibe et al.ll2012l for the rig- 
orous expression) with po ~ p{Ro,0), and s > 0. 

CsiR) 



H{R) 



n{R) 



(10) 



2.3 The simulations 

The simulations we have performed are listed in Table [T] 
The first column gives the name of the simulation, the sec- 
ond and third columns display the values of the two artificial 
viscosity (AV) parameters a and /3, the fourth the number 
A'" of SPH particles used for sampling the disc, the fifth the 
kind of artificial viscosity. Simulations can be divided in five 
sets. Each simulation can belong to more than one set. Sets 
are displayed in the last column of the table. Simulations in 
Set A and Set B allow to study the effect o f changing the ar- 
tificial viscosity parameters for the lMG83l and for the lMu96l 
artificial viscosity respectively. With simulations in Set C it 
is possible to study the effect of the different artificial vis- 
cosity implementations. In Set D are collected simulations 
designed for studying the effect of resolution. Finally, sim- 
ulations in Set E are used for testing the Ocont-a relations 
presented in Eqs.[7]and|8l The role of each set is summarised 
in Table [g] 

Midplane and vertical cross sections of the density pro- 
files of four of the performed simulations are shown in Fig. [l] 
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Table 1. Simulation list 



SI (0,1,0) 



Name 


a 


li 


N 




AV 


Set 


SI 


0.1 


0.0 


2- 


105 


MG83 


A,C 


S2 


0.1 


0.2 


2- 


105 


MG83 


A 


S3 


0.1 


0.5 


2- 


105 


MG83 


A 


S4 


0.1 


2.0 


2- 


105 


MG83 


A 


S5 


0.1 


10.0 


2- 


105 


MG83 


A 


S6 


1.0 


0.0 


2- 


105 


MG83 


A,C 


S7 


1.0 


2.0 


2- 


105 


MG83 


A,D 


S8 


1.0 


5.0 


2- 


105 


MG83 


A,D 


S9 


1.0 


10.0 


2- 


105 


MG83 


A 


SIO 


2.0 


0.0 


2- 


105 


MG83 


A,E 


Sll 


2.0 


4.0 


2- 


105 


MG83 


A 


S12 


2.0 


10.0 


2- 


105 


MG83 


A 


S13 


5.0 


0.0 


2- 


105 


MG83 


A 


S14 


5.0 


2.0 


2- 


105 


MG83 


A 


S15 


5.0 


10.0 


2- 


105 


MG83 


A 


S16 


0.1 


0.0 


2- 


105 


Mu96 


B,C 


S17 


1.0 


0.0 


2- 


105 


Mu96 


B,C,D 


S18 


2.0 


0.0 


2- 


105 


Mu96 


B 


S19 


1.0 


2.0 


5' 


lO'i 


MG83 


D 


S20 


1.0 


2.0 


1 ■ 


10*^ 


MG83 


D,C 


S21 


1.0 


5.0 


5' 


104 


MG83 


D 


S22 


1.0 


5.0 


1 ■ 


105 


MG83 


D 


S23 


1.0 


5.0 


5' 


105 


MG83 


D 


S24 


1.0 


5.0 


1 ■ 


10^ 


MG83 


D 


S25 


1.0 


0.0 


1 ■ 


10*^ 


Mu96 


D,C,E 


S26 


5.0 


0.0 


1 ■ 


lO'' 


Mu96 


E 


S27 


5.0 


0.0 


2- 


105 


Mu96 


E 



mm 




S2 (0.1,0.2) 



S6 (1,0) 



Figure 1. Disc morphology: volume gas density in midplane and 
vertical cross sections for simulations SI, S2, S6, S8 (from top 
to bottom and from left to right). The values of the two AV 
parameters [a, j3) are displayed for each case. 



erage quantities. In this paper we perform both the standard 
azimuthal and vertical averages {■)ez and the azimuthal av- 
erages only {•)e, in order to characterise also the vertical 
extension. 



Table 2. Simulation sets 



Set Description 



A Changing a and /3 in^G83 artificial viscosity 

B Changing a in ^u96 artificial viscosity 

C Changing the artificial viscosity model (MG83, ^u96) 

D Changing A^ 

E Testing the Ocont-o relations (Eqs.[7]and[8ll 



3 DIAGNOSTICS AND REFERENCE VALUES 
FOR PPD 

We consider two kinds of diagnostics which refer to the 
global flow and to fluctuations of several quantities, with 
particular focus on density and velocity field. We first define 
the diagnostic and then give the expected values in some 
known cases. In the following, the symbol (•) represents av- 



3.1 Diagnostic for the global disc flow 

Here, wc present the selected quantities used to characterise 
the global fiow. 

(i) Mass distribution. We look at: 

(a) the radial profile of the surface density E(i?), 

(b) the density distribution p{Rs,z) in the direction 
perpendicular to the disc midplane, for a selected radius 
Rs. 

(ii) Velocity structure. We focus on: 

(a) the radial velocity maps vr{R, z), where the radial 
component of the velocity field is azimuthally averaged. 

(b) the local Mach number of the flow in the midplane 
Ma(7?) = v{R)/cs{R), where v is the modulus of the local 
velocity field. 

(iii) Macroscopic turbulence signatures. Four items are 
analysed. 
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(a) The mass accretion rate onto the central star Mq: 
we are particularly interested in this quantity for two rea- 
sons: (1) it gives indirect information about the turbu- 
lent viscosity coefficient and (2) it is constrained by ob- 
servations. For PPD around T Tauri stars the measured 
values of mass ac cretion rates onto the central star are 
M » 10"* M0 (JHartmann et all Il998l : [Andrews et al.1 
l2009r). they can be reproduced by a v alue ass ~ 10 ~^ 

(JHartmann et al.lll998l : iKing et al.ll2007h . 

(b) Since Mo gives an information restricted to the very 
inner region of the disc (where our inner boundary condi- 
tion is free) , we also analyse the radial profile of the local 
mass accretion rate: M{R) — — 27r7?E{w_R)ez, averaged in 
the azimuthal and vertical direction. 

(c) The effective viscosity Ocff of the gaseous disc (we 
express the visc osity in terms of the dimensionless param- 
eter Qss of the |SS73| parametrization) that characterises 
the simulated disc. 

The effective viscosity in the simulated disc model is 
estimated by means of fits to the two-dimensional analytic 
models of accretion discs (see details in Appendix |B]). We 
use 2D analytic models because they are well suited for the 
azimuthal symmetry present in our simulated discs. For 
this reason we call Q2d this way of estimating the effective 
viscosity a^B. In particular, we derive a2D by comparing 
the vertical profile of the radial velocity of the gas flow 
to the expression derived by analytic Osg-disc models for 
different radial positions. We r ewrite the vertical profile 
of th e radial velocity, given by ITL02| and iFromang et al.l 
I2OIII . highlighting the scale height of the disc H{R) at 
radial position R: 



Vr 



Q2D H{R) 
2 Rq 



(15) 



6p + g - 3 - (9 - 5g) 



H{R) 



iLodato fc Pried (|2010|) estimated Qch by fit of the sur- 
face density profile of ID viscous accretion discs (jPringlg 
|l98l|) and found good agreement with Ocont.Muge in SPH 
simulations of warped discs when the necessa ry equ iva- 
lence conditions (large number of particles and IMu96| ar- 
tificial viscosity, see Sect. 12.131) ^"^^ met. 

(d) Once we know the effective viscosity of the disc, we 
can derive the corresponding Reynolds number: 



Ree 



Ma 

Oeff 



(16) 



3.2 Diagnostic for fluctuations 

For each simulated disc we have studied both the magnitude 
and the structure of the fluctuations present in the disc with 
particular attention to the velocity field. Since current ob- 
servations have not yet reached the resolution necessary to 
directly detect turbulence and study its features in proto- 
pl anetary discs (some pr ogress has been recent ly claimed 
bv lHughes et al.ll201ll and iGuilloteau et al.ll2012l . who mea- 
sured some turbulent velocities in mm observations) the re- 
sulting properties of density and velocity fiuctuations have 
been compared to the typical behaviour observed in turbu- 



lence experiments and to results from grid- and particle- 
based simulations available in the literature. 

The components of the velocity field are identified by Vi 
and those of the fluctuating velocity field are ut = Vi — (vi) 
with i = R,0, z. 



3.2.1 Magnitude of velocity fluctuations 

Concerning the magnitude of velocity fiuctuations, we focus 
on the turbulent viscosity coefficient (point Ml), the diffu- 
sion coefficients (point M2) and the Mach number of velocity 
fluctuations (point M3). 

(Ml) The turbulent viscosity coefficient vt is often de- 
rived from Reynolds averages of Navier-S tokes equa tions 
with the turbulent viscosity hypothesis (e.g. |Popell2000r ) and 
is thus related to Reynolds stresses: 

a{{ve)e/Ry 



i^T = —{uRUe)e I 



R- 



dR 



(17) 



(see Appendix lAl[) . 

Here we call ors the corresponding [SS73 parameter for a 
disc in quasi-Keplerian rotation, it is given by: 

2 {uRUe)e 

3 c2 ' 



aB.s{R,z) 



(18) 



where Eq. [17] has been combined with the |SS73| relation 
vt = assCsH (see Appendix I A2[) . For each radial and ver- 
tical position, averages are performed spatially along the 
azimuthal direction in order to be able to derive the radial 
and vertical dependence of qrs. 

The corresponding turbulent Reynolds number is: 

Ma 



Rbt — 



ans 



(19) 



In MHD simulations of accretion discs, Oss ~ 10 ^ 
and includes both the Reynolds and Maxw ell stresses (e.g. 
iFromang fc NelsonI [2OO9I : iFlock et alll2012l) . in simulations 
of self-gravitating discs, Oss ~ 10~^ and i ncludes both the 
Reyn olds and gravitational stresses (e.g. iRice et al.l |2005| . 
[20i3). 

(M2) Turbulent flows are often approximated by means 
of models where turbulence is described as a diffusive pro- 
cess. Here we want to determine if diffusion is present 
in our simulations. To this end the method used by 
IFromang fc Papaloizoul (|2006l ) in order to calculate the tur- 
bulent diffusion coefficient Dt is well suited for Lagrangian 
codes such as SPH: 

DT{t) = / SMt')dt' (20) 

Jo 

where Szz is the vertical velocity correlation function: 



S,,{t) = {v4z,t)v4zo,0)}, 



(21) 



with zq representing the vertical position of a given particle 
at time i = and z the vertical position of the same particle 
at time t; in a similar way Vz{zo,0) is the velocity of the 
same particle at time t = and Vz{z, t) at time t. In order to 
compute the diffusion coefficient at a selected radial location 
Rb, we average the diffusion coefficient of single particles 
that were at Rs at the beginning of the simulation. 

A diffusion coefficien t Dt/ (cbH) « 5 ■ lO"'^ is found by 
IFromang fc Papaloizoul (|2006r ) in MHD local shearing box 
stratified simulations. 
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(M3) The strength of velocity fluctuations and therefore 
their subsonic or supersonic behaviour can be determined 
by the corresponding Mach number, defined as the ratio of 
the modulus of velocity fluctuations to the gas sound speed: 



Maf^^. 

Ca 



(22) 



We have computed azimuthal-only and azimuthal and verti- 
cal averages in order to present R-z maps and radial proflles. 
In the global MHD simul ations of stratifle d discs of 
iFromang fc NelsonI (|2006l ') and lFlock <^~^ (|201lf ). the Mach 
number of velocity fluctuations is Maf « 0.1 around the mid- 
plane and Maf ~ 0.4 in the surface layers. 



3.2.2 Structure of velocity fluctuations 

Concerning the structure of velocity fluctuations, we focus 
on the anisotropy vector (point SI), on the power spectrum 
(point S2) and on the possible presence of intermittency 
(point S3). 

(SI) In order to determine the isotropic or anisotropic 
character of velocity fluctuations, we deflne the anisotropy 
vector 5 = [5r, 50,5z) in terms of the velocity dispersion in 
the three sp atial directio ns, in analogy with galactic dynam- 
ics (see e.g. lBertin|[200d '): 



„2 I „2 



(23) 



where the square of the velocity dispersion a^ — {{vi — 
(vi))'^) = {u\) = {v1) — (vi)^ and i, j and k refer to the 
three cylindrical components. Systems with anisotropy in 
the i direction are characterised by a velocity dispersion of 
the i component that is larger than that of the other two 
components: at 3> o-j and cri^ a^. Therefore, the i compo- 
nent of the anisotropy vector will he Si Ki 2. Averages are 
performed along both the vertical and the azimuthal direc- 
tion. 

From th e plo ts of iFromang fc NelsonI (|2006l ) and 
iFlock et al.l (|201ll ) we can deduce that or S> ere « az, 
therefore MHD discs are radially anisotropic and the cor- 
responding anisotropy vector is 5 ~ (2,0,0). In contrast, 
in experiments of turbulence produced by round jets a cor- 
responding azimuthal anisotropy 5 ~ (0, 2, 0) has been ob- 
served (e.g. lPopellioOol ). 

(S2) The power spectrum is computed along a ring of se- 
lected radius Rs, centred on the star and located at the 
selected vertical position Zs (see Appendix I A3|) : 



P,{R„Zs;k) = \Mk)\" 



(24) 



with Mi(fc) the Fourier transform of the component i of the 
velocity fluctuation vector and k = 1/A the wavenumber 
corresponding to length scale A. We analyse the power spec- 
trum of the velocity fleld in order to determine if an energy 
cascade, characteristic of turbulent systems, is present and 
if differences related to the position inside the disc are rele- 
vant. 

(S3) Hi ghly turbulen t flows are characterised by intermit- 
tency fe.g. lFrischlll996r ). This feature implies non-Gaussian 
probability distribution functions (PDFs) with correspond- 
ing non-Gaussian higher order moments. Here we focus on 



the density and acceleration PDFs and on their 3"^"^ {skew- 
ness, noted S) and 4*'' {kurtosis, noted K) order moments 
(see Appendix I A3fl . 



Data are not available for the power spectrum of fluc- 
tuations and for the phenomenon of intermittency in accre- 
tion discs. In the classical Kolmogorov theory of incompress- 
ible homogeneous turbulence the power spectrum of velocity 
fluctuations has a power law form P{k) ~ k" with s = —5/3. 
Recent simulation s of supersonic c ompressible gas found a 
slope of s « —2 (jPrice fc Federrath 201Qi , and references 
therein). Intermittency, highlighted by non-Gaussian PDFs 
and higher order moments, is experimentally found in very 
high Reynolds number incompressible flows. In SPH simu- 
lations of a s imple weakly compressible periodic shear flow, 
lEllero et af] (|2010l ') flnd that the PDF of the particle accel- 
eration is in good agreement with non-Gaussian statistics 
observed experimentally for incompressible flows. 

In s imulations of supersonic compressible turbulence 
(see e.g. IPrice fc FederrathluOlOl ). log-normal distributions 
of the density PDF have been observed: 



p{x) = 



^7, 



: exp 



{x-{x)r 



2al 



(25) 



with X — ln(p/{p)) and a width controlled by the Mach 
number of the fluctuations 



al = In (1 + &^Ma^) 



(26) 



where b ~ 0.5. The expressions of the skewness and kurtosis 
in terms of ap for the log-normal distribution are 



S 



(e^p +2] Ve 



K = e^^p + 2e'^''p + 3e^''p 



(27) 



When (jp — >■ the log-normal distribution tends to a Gaus- 
sian distribution and the skewness and kurtosis tend to their 
gaussian values (0 and 3 respectively). 

The results of the application of these diagnostics to the 
simulations we performed are described in Sects 2] E] and [B] 



4 THE DISC GLOBAL FLOW 

We now consider how the main properties of the global disc 
flow (mass distribution, velocity structure, mass accretion 
rate and effective disc viscosity) depend on the strength of 
artiflcial viscosity. The convergence of each quantity is also 
studied. 



4.1 Mass distribution 

The position of the peak of the surface density proflle E(_R) 
moves towards larger radial distances and its maximum de- 
creases (see left panel of Fig. [5}, when the two AV parame- 
ters a and /3 are increased: such behaviour is in agreement 
with the viscous spreading mechanism that is expected here 
due to the existing correlation between artiflcial and physi- 
cal viscosity (see Sect. l2.L2|) . 

An increase in resolution (see right panel of Fig.[2|) has a 
qualitatively similar effect to a decrease in viscosity: the sur- 
face density peak moves towards the centre of the disc, be- 
coming higher (note also that the curve becomes less noisy). 
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Figure 2. Mass distribution. Dependence of the surface density 
S(iJ) and vertical density profile p{Rs = ^,z) on the AV param- 
eters (left: simulations S6, S7, S8, S9 with N = 2 ■ W^ particles) 
and on the resolution (right: simulations S21, S8, S23, S24 with 
(a,/3)=(l,5)). Values are given in code units. 



In the central and external regions the projected density pro- 
file converges already for a number of particles as small as 
5- 10*. In the inner region [R < 1Q~^ Rout) , a trend towards 
convergence is present, but a higher number of particles, 
N > 10®, is necessary. 

The vertical density profile p{Rs, z) at a selected radius 
Rb (see insets in Fig. [2] for Rs ~ 1) follows the expected 
Gaussian distribution of Eq. [9] showing that the desired 
profile is well reproduced by the sampling procedure and 
by numerical thermalisation. It is not significantly affected 
by changing a, 13 or N. 



4.2 Velocity structure 

4-. 2.1 Radial velocity maps 

Maps of the azimuthally averaged radial velocity {vR)g in 
the R-z plane are displayed in Figs. [3] and |4] for changing 
artificial viscosity and resolution respectively. 

Effect of artificial viscosity parameters. A transition 
from a chaotic to a regular accretion pattern is clearly ob- 
served in Fig. [3] when a is increased for a given /3 (panels 
from top to bottom) and as well as when /3 is increased for 
a given a (panels from left to right). 

At large radii, a well extended decretion region domi- 
nates in all models, as expected from the conservation of an- 
gular momentum. At smaller radii, instead, the structure of 
the gas flow strongly depends on the values of the AV param- 
eters. A decretion flow located around the midplane, asso- 
ciated with an accretion flow in thin surface layers, appears 
when (a, /3) are increased. Such accretion layers become 
thicker for larger AV parameters. Therefore, the gas flow is 
characterised by the phenomenon of meridional circulation. 
This mechanism has bee n foun d in some two-dimensional 
viscous disc mo dels ( e.g. ITL02| that extends the standard 
one dimensional I SS 731 models to two dimensions: radial and 
vertical) but is not reproduced by MHD d isc simulations 
where turbulence i s induced by the MRI (|Fromang et al.l 
I2OIII : |jacauetll2013l ) . The presence of meridional circulation 
has important implications for the topic of particle mixing 
in PPD and is one of the possible mechanisms able to ex- 
plain the presence of crystalline solid particles in the outer 
regions of T Tauri stars recently observed by Spitzer (e.g. 
ISargent et al.ll2009 ). Finally, in the very inner region of the 



disc a small and well defined accretion region is present in 
all models. 

Effect of artificial viscosity implementation. The two 
panels in the bottom row of Fig. [4] show the structure of the 
radial velocity flow of two discs at different resolution sim- 
ulated with the lMuQa artificial viscosity. The resulting flow 
structure can be directly compared with the panel above, 
which represents the flow structure in an e quivalen t disc at 
the same resolution but simulated with the|MG83J artificial 
viscosity. In the simulations at lower resolution (N = 2- 10'''), 
the meridi onal circulation pattern is better reproduced by 
the IMG83I implementation. However, in the more resolved 
discs {N = 10*^) the two patterns are very s imilar. The two 
small accretion regions present in the outer iMuQa discs are 
due to a longer relaxation time-scale for this particular ar- 
tificial viscosity implementation (in longer evolution simu- 
lations, not shown here, we have ob served a net outward 
gas fiow in agreement with the lMG83l case and the expected 
viscous spreading). 

Effect of resolution . The first row in Fig. [4] shows how 
the radial velocity structure has already converged after N ~ 
5 • 1 0^ par ticles. 

ITL02| found that the outflow zone around the midplane 
shrinks as the radial volume density gradient is reduced (s — >■ 
0) and when the condition p + q < 2, that guarantees a net 
accretion of the gas onto the star, is violated. Therefore, we 
expect that the thickness of the accretion layers depends 
on the values {p, q) of the surface density and temperature 
radial profiles of the gaseous disc model. However, a detailed 
study of the dependence of the structure of the accretion 
layers on the gas density and temperature profiles and the 
correlation with the central mass accretion rate is out of the 
scope of this paper and will be addressed in a future work. 
For the initial profile of the disc model used in this work 
we have p + q = 2.25, which implies net decretion but is 
very close to the limit between net accretion and decretion. 
The surface density of the simulated discs after numerical 
relaxation reaches a smooth profile that, in the inner region, 
is slightly shallower (implying a slightly smaller p) than the 
initially imposed power law, which is characterised by an 
inner sharp edge (see Sect. 12. 2|) . This moves the model in 
the net accretion region and explains the presence of the 
observed thin accretion layers. 



4-2.2 Mach number of the flow 

The gas in the disc is supersonic. In fact, since the gas 
moves on quasi-Keplerian orbits, the Mach number is ap- 
proximated by Ma ~ Uk/cs. With the adopted sound speed 
profile it becomes: Ma « R^''~^''^ /cbo. For the values used 
here Ma ~ 20-R~^'*, which is approximately in the range of 
values 15-25. 



4.3 Mass accretion rate 

The mass accreted onto the central star (Fig. [S| and the 
radial profiles of the mass accretion rate (Figs. [6] and [7|) are 
two macroscopic signatures of turbulence. 
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Figure 3. Radial velocity structure of the gaseous disc: changing 
AV parameters. Maps of the azimuthally averaged radial velocity 
{vji {R, z))g in th e meridi an plane for simulations with N = 2- 10^ 
particles and the JMGSSl artificial viscosity implementation. The 
colorbar gives the radial velocity; negative values correspond to 
gas inflow and positive values to outflow. Values of (a,/3) arc given 
in each subplot. 



: — i i r 

(1 ,5) - M.5S4 






I ° 

'^■0.4 

-0.8 



-t-+ 



(l,0)Mu96-N=2o 
S 0.4 - 



-<!:■: 



'*'. 



-< 



<,'J*^ 



0.8 
5 0.4 

< 

g 
N -0.4 

-0.8 



If"'. 



-I— I- 



H 0.001 


0.001 



■ 0-001 



<-Ti 



H 



12 3 
R(IOOAU) 



1 2 3 
R(100 AU) 



Figure 4. Radial velocity structure of the gaseous disc: chang- 
ing resolution and artificial viscosity implementation. The lay- 
out is the same as in Fig, \3\ The top row shows how the radial 
vel ocity st ructure changes with resolution for simulations with 
the lMGSSi artificial viscosity, the specific case of simulations with 
(a,/J)=(l,5) has been chosen as an example: panels from left to 
right correspond respectively to A'"=5 ■ 10*, 2 ■ 10^ , 5- 10^ , 1 • lO*". 
The bottom row refers to two simulations with the|Mu96| artificial 
viscosity implementation (with a = 1) and with N=2 ■ 10^ (left 
panel) and N=l - 10® (right panel). 



4-. 3.1 Central mass accretion rate 

Effect of artificial viscosity. The left panel in Fig, [5] shows 
that Mo inc reases with larger a and/or larger /3 and in the 
case of the IMu96| artificial viscosity implementation. The 
effect of Q is larger, as expected since it controls the first 
order in the artificial viscosity term. 

Effect of resolution . The right panel in Fig. [5] shows 
that Mo decreases with increasing resolution. Such an effect 
is due to a decrease of particle noise at higher resolution, 
which decreases the numerical dissipation. This is linked to 
the dependence of the effective viscosity on resolution (see 
Sect.OJ. In fact, as will be shown in Sect. 16.3.11 the PDFs 
of several physical quantities are narrower at higher resolu- 
tion, due to the reduced noise. Concerning the central mass 
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Figure 5. Central mass accretion rate Mq. Left panel: depen- 
dence on the AV parameters for a given number of particles 
(N = 2 ■ 10^), Right panel: dependence on resolution. The mass 
accretion rate is displayed in units of Mq yr^^. 
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Figure 6. Gas flow in the disc. Surface density radial profile (top 
panel) and mass accretion rate radial profile (bottom panel) for 
simulation S8, 



accretion rate we look at the part of the radial velocity distri- 
bution corresponding to negative velocities, since only fluid 
elements with negative radial velocity are responsible for ac- 
cretion. The radial velocity distribution at higher resolution 
presents a smaller spread than in the lower resolution case. 
In addition, SPH particles in low resolution simulations have 
a larger mass than in higher resolution simulations. There- 
fore, at lower resolution there is a larger number of more 
massive SPH particles with higher inward radial velocity, im- 
plying higher accretion rates. The trend of the curve shows 
that the central accretion rate is converging, however con- 
vergence is not yet reached at one million particles. 

For all the considered combinations of {a,/3) and A'', we 
find mass accretion rates consistent with the values observed 
for PPD. In addition, we observe a correlation between the 
trend of Afo with respect to (a, /3) and the increase of thick- 
ness of the accreting layers (net accreting mass flux) ob- 
served in Fig. 13] Similarly, the trend with respect to A^ cor- 
relates with the decrease of thickness of the accreting layers 
at higher resolution, observed in Fig. |4] 
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Figure 7. Radial profiles of the mass accretion rate. The effect of 
artificial viscosity is shown for a = 0.1 (top panel) and 1 (middle 
panel) and changing /3, for simulations with N = 2 ■ 10^ parti- 
cles. The effect of resolution is displayed in the bottom panel for 
simulations with (a, /3) = (1,5). 



4-3.2 Radial profile of the mass accretion rate 

The central mass accretion rate Mq gives information only 
about the region of the disc that is close to the star. In or- 
der to have a global picture of the mass flux in the disc, it is 
interesting to look at the radial profile of the vertically inte- 
grated mass accretion rate M[R). As an example, the profile 
observed in simulation S8 is displayed in the bottom panel of 
Fig. [S] The flow is characterised by an inner accreting region 
separated, at a radius we call Rt, from an outer decreting 
region with small rates at intermediate radii and larger rates 
at large radii. It can be understood by a comparison to the 
surface density profile (top panel of Fig. O . In the organized 
fiow of model S8, Rt corresponds to the radial location of 
the surface density peak, which moves outwards with time 
due to viscous spreading, as the density maximum decreases. 
The density peak behaves as a reservoir of gas that supplies 
mass both to the inwards and outwards flows, resulting in 
an apparent null mass accretion rate at Rt . The varying rate 
in the decretion region is due to the vertical structure of the 
flow: Fig. 13] shows accretion throughout the disc height for 
the inner disc, accretion in the midplane with decretion in 
the surface layers at intermediate radii resulting in a slightly 
negative radial mass flux, and decretion throughout the disc 
height in the very outer regions. 

Effect of artiflcial viscosity. For low values of a and /3 
fluctuations are so large that the identiflcation of an orga- 
nized flow of the fluid is not possible. As shown by the top 
panel of Fig. [T] the mass accretion rate proflle is very noisy 
with a narrow peak in the inner disc, an average value close 
to zero in the intermediate regions and negative in the outer 
disc. The peak of the surface density proflle decreases with 
time but less than in the higher {a, /3) case and it does not 
change position with time. In the vr map (see the top left 
corner of Fig. [3} negative and positive regions are clearly 
mixed. 



Increasing the two AV parameters, the proflle becomes 
more regular and smooth, following the structure described 
above for the speciflc simulation S8. For a = 1 the proflle 
tends t o conv erge when /3 is increased (see the middle panel) . 
In the iMuQa case we observe similar radial proflles for the 
mass accretion rate (not shown). 

Effect of resolution . As shown in the bottom panel of 
Fig.[7|for the simulations with (a,/3) = (1, 5) and increasing 
resolution, the peak of the radial proflle of the mass accre- 
tion rate moves inwards and converges for N > h ■ 10''', in 
agreement with the convergence observed for the radial ve- 
locity maps. The reason of the decrease of the peak in higher 
resolution simulations is the same that explains the similar 
decreases observed for the central mass accretion rate (see 
Sect. HXTj) . 



4.4 The effective viscosity of the SPH disc: Q2d 

In order to estimate the effective viscosity of the gaseous 
disc, we determine the 0120 parameter, which is derived by 
fltting the vertical proflle of the radial velocity expected for a 
viscous axisymmetric accretion disc, given by Eg. 1151 to that 
derived from the simulations (see Appendix [B] for the flt- 
ting procedure). Two examples, characterised by adequately 
large AV parameters, are shown in Fig. [H] where the verti- 
cal profile of the radial velocity at radial position _R = 1 is 
displayed for simulation S7 (top panel) and for simulation 
S20 (bottom panel). 

In Fig. [9] we present the values of a2D obtained by av- 
eraging the result of the fit at R=l and R=2 for the set of 
simulated discs. The values are averaged in the radial range 
[1,2]. In simulations with a = 0.1 the vertical profile of ra- 
dial velocity cannot be described by the profile in Eq. [TJ] 
because the high noise level makes the fit impossible. We 
find that simulation profiles can be fitted by a value of 020 
which is roughly the same for the different radial locations 
we considered. This shows that the str ucture of SPH mod- 
els is consistent with the 2D version of |SS73| accretion disc 
models, characterised by a constant a^s value. 

Effect of artificial viscosity. In the range of AV parame- 
ters were the determinations of the effective viscosity is pos- 
sible (a > 0.1), the scaling of 020 with the AV parameters 
(increase with larger a and/or /3 in the iMGSa implementa- 
tion and increase with larger a for the IMuQGI implementa- 
tion) has the same qualitative behaviour as for the central 
mass accretion rate (see Fig. [5| . Such a match is what we 
expect, since a more viscous disc is characterised by a larger 
central accretion rate. The weaker dependence on j3 is in 
agreement with the fact that B is the parameter controlling 
the second order term of the iMGSa artificial viscosity (see 
Eq.SJ. 

Effect of resolution . The trend of the curves in the right 
panel of Fig. [9] is qualitatively similar to that of the central 
mass accretion rate (see Fig. [SJ and shows that the effec- 
tive viscosity is converging. However one million particles is 
still not enough in order to reach convergence. The slightly 
higher resolution required for the convergence of Qf2D and Mo 
with respect to that required by the radial velocity maps and 
by the mass accretion rate radial profile (A*' ~ 5 • 10^) is due 
to the fact that the former two parameters are calculated 
from a smaller subset of simulation particles, the numerical 
noise is therefore larger. 
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Figure 8. Determination of a2T>- The vertical profile of the radial 
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simulation S7, bottom; simulation S20. 
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Figure 9. Results for a2D- The parameter has been determined 
by averaging the values fitted at R=l and 2 and only for simula- 
tions where the noise is low enough to allow the fit. 




Figure 10. The effective Reynolds number: Re2i3 . Derived from 
the effective viscosity presented in Fig. [9] 



to Q, /3 and A'', described in the last paragraph, we found 
larger Re numbers in less viscous and more resolved discs. 

Simulations with one million particles and the JMGSSl 
artificial viscosity are closer to, but sill lower than the value 
Reofi « 3 • 10^, which is considered as a limit for the onset 
of t he turbulen t regime in a box with periodic boundaries 
(see |Pricell2012l ). Therefore, the dynamics of the discs pre- 
sented in this paper is not expected to be dominated by 
fully-developed turbulence. 



4.4-2 The a2D-Qcont connection 

Figure [TT1 shows the relation between the effective viscosity 
020 (computed in the simulations) and the Ocont viscosity 
(expected when the continuum limit of the SPH equations is 
taken) for eight of the performed simulations with P — (S6, 
SIO, S13, S17, S18, S25, S27). Simulations SI and S16 (with 
a = 0.1) are too noisy to allow the determination of 020- 
Since the ratio (h) /H in the expression of Qcont depends on 
the radial location inside the disc (see Eqs.[7]and[8]), we have 
evaluated locally at the radial position R — 1 both acont and 
a2D, whose values are given in Tabled 

We find that both th e point s concerning the lMu96l case 
and those concerning the IMG83| case follow the appropriate 
analytic relation. We therefore confirm the necessity of re- 
ducing the coefficient in the Qcont relation i n Eq. |8] for the 



lMG83l case as shown b y I Mem fc Bate! (|2012l ) and motivated 
by the fact that in the JMGSSl implementation the artificial 
viscosity is applied only to approaching particles. Finally, no 
significant difference is observed when resolution is changed 
from 2 • 10^ to 10^ particles, suggesting that already with 
2 • 10^ particles the continuum approximation can be con- 
sidered valid. 



The result of this subsection is that the effective viscos- 
ity 020 in all the SPH models considered here equals a few 
10~^, in agreement with the order of magnitude deduced 
from observations. 



4-4-^ Reynolds number 

Once the effective viscosity of the disc is determined 
(ocff ~ 020), it becomes possible to derive the correspond- 
ing Reynolds number Ke^s ~ Re2D = Ma/a2D, shown in 
Fig. 1101 In accordance with the trend of Q2d with respect 



5 THE MAGNITUDE OF SPH 
FLUCTUATIONS 

5.1 The Reynolds stress contribution to the SPH 
disc viscosity: a-Rs 

We now proceed with the determination of a-RS , which is the 
ass parameter derived from the Reynolds stress, in analogy 
with turbulence studies and with MHD or self-gravitating 
disc simulations. Radial and azimuthal fluctuations are com- 
puted from the average velocity field determined locally us- 
ing a R-z grid. 

The radial profile of the vertically averaged aas is 
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Figure 11. Comparison of the the effective viscosity of the disc 
02D with the co rresponding Qcont formula: simulations S6, SIO 
and S13 | |MG83|. N=2 - 10^) S17, S18 and S27 l lMuQd N=2 ■ 
10^) and S25 and S26 JMuQd Af=10'^). The displayed values are 
computed at the radial position R = I. The dotted line represents 
identity between Ocont and «2D- Error bars are comparable to 
the symbol sizes and are not displayed, their values are given in 
Table [3l 



Table 3. Error bars associated to 0.20 produced by the fitting 
procedure. 



Simulation 



^cont 



S6 


0.0301 


0.0490 ± 0.0063 


SIO 


0.0597 


0.0668 ± 0.0040 


S13 


0.1476 


0.1453 ± 0.0023 


S17 


0.0601 


0.0759 ± 0.0043 


S18 


0.1204 


0.1249 ± 0.0033 


S27 


0.3009 


0.2949 ± 0.0075 


S25 


0.0372 


0.0500 ± 0.0013 


S26 


0.1852 


0.1776 ± 0.0009 



shown in the left panel of Fig. 1121 where the effect of chang- 
ing a and /? is highlighted, and of Fig. 1131 where the effect 
of resolution is shown. In all simulations, the coefficient qrs 
is characterised by a similar profile: it tends to be approxi- 
mately constant throughout the radial extension of the disc, 
with an increase both in the central and in the external ra- 
dial region. 




Figure 12. Radial (left) and vertical (right) profiles o/qrs- Ef- 
fect of changing the AV parameters. 
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Figure 13. Radial (left) and vertical (right) profiles of ajt^s- Ef- 
fect of changing the resolution. 
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Figure 14. Vertical profiles o/a^g. Comparison at three differ- 
ent radial location in the disc for simulation S24. 



The right panel of both figures shows the radially av- 
eraged vertical profile (as in Fig. [8l the semi-thickness 
of the disc is computed from H{R)=Ho{R/Ro)^^~''^^^ , 
see Sect. I2.2|l : ors is approximately constant around the 
midplane and then increases with height. This behaviour 
qualitatively agre es with that observed in MHD discs by 
iFIock et al.l (|201ll ). even if we observe a broader minimum. 
In addition, we have found that the vertical profile is ap- 
proximately the same at different radial positions, as shown 
in Fig. [T4]for simulation S24. 

Effect of artificial viscosity. The radial and vertical pro- 



files of aRs are very sensitive to the AV parameters a and (3. 
In fact there is a sharp change in the average value, which 
decreases from positive to negative values, when the two AV 
parameters are increased. This transition is clearly visible 
in both panels of Fig. [12] for a — 1: the increase of /3 from 
to 2 shifts the profile of qrs towards smaller values, moving 
it in the negative region. A further increase of /3 leads to 
smaller changes, but the behaviour inverts: higher artificial 
viscosity now tends to increase the value of the coefficient 
Qrs (which corresponds to a decrease of its absolute value, 
see left panel) and to extend the constant region around the 
midplane (see right panel). We find the same trend for other 
values of a (not shown). 

The change of sign of ors is due to the mechanism of 
meridional circulation that is correctly resolved when ade- 
quately high AV parameters are used. 

Effect of resolution . For a given combination of AV pa- 
rameters, the radial and vertical profiles of the q:rs coeffi- 
cient do not change with resolution. The only effect of in- 
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Figure 15. Time evolution of the diffusion coefficient. D'y/{csH) 
computed at ij = 2. Left: low artificial viscosity simulations (SI, 
S2, S3 , S4). Right: higher artificia l visco sity (S7, S8, S24 with 
[MG83 implementation and S17with[Mu9g implementation). Note 
the different vertical scales. 
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Figure 16. Results for Dt- Left: effect of artificial viscosity, 
right: effect of resolution. 



creasing the number of particles of the simulation is a re- 
duction of the noise, which implies smoother profiles. 

We note that the computation of Qrs is very sensitive 
to the correct determination of the average velocity field 
of the flow, therefore only in simulations with adequately 
high values of a and j3, where numerical noise and particle 
disorder are lower, can it be trusted. The absolute value of 
aas for the present simulations is |qrs| ~ 10"'^. This is 
comparable to standard values found in MHD discs (see e.g. 
iFromang fc NelsonluOOSl ). The difference is that meridional 
circulation is present here and the negative value indicates 
decretion on average in the midplane, in contrast to the 
positive value of MHD discs (where meridional circulation 
is not seen), which is expected for accretion in the midplane. 

In conclusion, more viscous discs are characterised by 
a smaller modulus of cvrs and, for all the simulations we 
performed, we find the ratio of the Reynolds stress to the 
effective viscosity of the gaseous disc to be aRs/Q2D ~ 10^^. 
This means that velocity fluctuations are not the main com- 
ponent of the effective viscosity of the disc, in agreement 
with the low Reynolds number of the global flow in in the 
present simulations (see Sect. 14. 0]) . However, the turbulent 
Reynolds number, which is associated to the Reynolds stress, 
amounts to ReT ~ 10'* (since Ma ~ 10 and OfRs ~ 10^'^), 
close to the turbulent limit. Therefore, physical fluctuations, 
quantified by the small contribution of ors to the effective 
viscosity, can be expected to exhibit turbulent signatures. 
We will show in Sects. [5^ and [6^ that it is indeed the case. 

An important point to be highlighted is that the de- 
termination of ors is not a sufficient condition in order to 
determine the mass accretion rate of the disc, since other 
sources can be present in the model. 



5.2 Diffusion coefficient 

For each simulation, we have computed the diffusion coef- 
ficient from Eq. [20] at several radial positions in the disc. 
In Fig. [TJ] we show the evolution with time of the diffusion 
coefficient Dt, in units of CsH, measured at the intermedi- 
ate radial position R — 2 for some simulations belonging to 
Set A. Except for simulations from SI to S3, characterised 
by larger fluctuations, the diffusion coefficient regularly in- 
creases with time and converges to a well-defined value. 

For the standard combinations of parameters {a, /3) = 
(1,2), D t/{cbH) ~ 5 ■ 10""^ is in agree ment with the value 
found bv lFromang fc Papaloizoul (|2006l ) in their MHD local 
shearing box stratified simulations. The asymptotic values 



of the diffusion coefficient for each performed simulation are 
displayed in Fig. 1161 

Effect of artificial viscosity. The left panel of Fig. [16] 
shows that the coefficient decreases with increasing a and 
P (except simulations from SI to S3 which present a diffu- 
sion coefficient without a clear trend, due to the higher noise 
level) and in the 1Mu96| case (the artificial viscosity switch 
on approaching/receding particles is off, resulting in more 
viscous discs). This trend is due to a reduction of the noise 
that allows the onset of meridional circulation, which forces 
the fluid to be more localised in the vertical direction. 

Effect of resolution . From the right panel of Fig |16l we 
observe that the diffusion coefficient decreases with increas- 
ing resolution. The reason of this behaviour is the same that 
explains the qualitatively similar trend we have observed 
in Sect. 14.3.11 for the central mass accretion rate. In this 
case, the PDF of the vertical velocity Vz is broader in low 
resolution simulations with respect to more resolved sim- 
ulations. Therefore, larger vertical velocity are possible in 
less resolved discs with resulting higher values (see Eq. I21|) 
of the diffusion coefficient. The curve shows a converging 
trend, however convergence is not yet reached at one mil- 
lion particles. As for the central mass accretion rate and for 
the effective viscosity fit, convergence requires more reso- 
lution because the diffusion coefficient is computed from a 
smaller subset of particles than in the case of the radial ve- 
locity maps and the mass accretion rate profiles, where an 
average on the azimuthal component is taken. 

We conclude that the diffusive mechanism is correctly 
represented by intermediate/large AV parameters (a,/3): 
more viscous discs are less diffusive, due to the onset of 
meridional circulation and the decrease of aRs . In addition, 
for a given (q, /3) combination, less resolved discs are more 
diffusive: such numerical dependence is of the same type as 
the dependence of the central mass accretion rate and of 
the effective viscosity on the resolution of the simulation 
(see Sect. [1331) 

5.3 Mach number of velocity fluctuations 

The azimuthally averaged Mach number of the modulus of 
velocity fiuctuations is displayed in the R-z plane in the left 
panel of Fig.[T7]for simulation S23, as an example. The right 
panel shows the corresponding radial profile for each velocity 
component, vertically and azimuthally averaged. All simu- 
lations show qualitatively similar structures and profiles. 
Effect of artificial viscosity. The average value of Mat 
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Figure 17. Mach number of velocity fluctuations. Left: azimuthal 
average in tiie R-z plane, right: radial profile of the vertieal and 
azimuthal average for each velocity eomponent, in simulation S23. 




Figure 18. Results for Maf . Effect of changing the AV parame- 
ters and the AV switch (left panel) and resolution (right panel). 



in the radial range R € [1, 2.5] is displayed in Fig. 1181 From 
the left panel of the figure, we see that, except for simula- 
tion SI, whose fluctuations are close to the sonic regime, all 
simulations show subsonic fluctuations with Mach number 
signiflcantly smaller than unity. We observe a decrease of the 
average Mach number with larger AV parameters. For rela- 
tively high a and /?, at each given radius the Mach number 
increases from a minimum value in the midplane (Maf ~ 0.1) 
to a maximum value in the surface layers (Maf ~ 0.2), as dis- 
played by the left panel of Fig. 1171 This behaviour is qualita- 
tivel y in agreement with what i s found in MHD s imulations 
(e.g. iFromang fc NelsonI [JOOgI : iFlock et al.1 120111 ), however 
we observe smaller values. A global maximum Maf « 0.4 
is reached in the inner region, where vh < and the two 
accreting layers meet. In addition, two layers with locally 
higher fluctuations with Maf ~ 0.2 are also present at the 
interface between the surface accreting layers and the mid- 
plane outwards flow. 

Effect of resolution . Similarly, as a consequence of noise 
reduction, higher resolution leads to a decrease both of the 
Mach number and of the extension of the surface layers 
where it reaches the local maximum for a given radial po- 
sition (with the same qualitative behaviour of the accreting 
layers in the vr maps, see Figs. [3] and |3| , as expected. We 
have found a convergence of the vertically and azimuthally 
averaged Mach number for a number of particles Af ~ 5 • 10'' 
(see right panel of Fig lTHj) . in agreement with the value found 
for radial velocity maps (which are connected to radial ve- 
locity fluctuations and therefore to the radial part of the 
Mach number). 

Fluctuations tends to be dominant in the regions of the 
disc where accreting and decreting flows meet: the cone- 




Figure 19. Results for 5, the anisotropy vector of velocity fluctu- 
ations. Effec t of ch anging the AV parameters and the AV switch 
IIMGSa and iMuQa ) , values are averaged vertically, azimuthally 
and in the radial domain [1, 2.5]. Note that values of the azimuthal 
and vertical anisotropy (middle and right panel) for /9 = become 
highly negative, they fall outside of the displayed area. 



like regions between the accreting surface layers and the 
decreting midplane layer, and the region around the inner 
disc edge. 



6 THE STRUCTURE OF SPH 
FLUCTUATIONS 

6.1 Anisotropy of velocity fluctuations 

We recall that when Sr = 2 the system is totally radially 
anisotropic, while for negative values the system is domi- 
nated by anisotropy in different directions. The same holds 
for the azimuthal and vertical components. The kind of 
anisotropy of the disc is displayed in Fig. [19] for simulations 
with changing a and (3. We observe that none of the simu- 
lated discs is totally anisotropic in one particular direction, 
but there is always a component that dominates the other 
two. 

Effect of artificial viscosity. Discs with smaller AV pa- 
rameters are dominated by radial anisotropy. For example, 
in the case of simulation SI with (a, P) — (0.1,0), we see from 
Fig. [T9]that 5r ~ 1.4, Se < —1 and 5^ < —1, it follows that 
5r is signiflcantly larger both than 5$ and than 5^ . When a 
and/or /3 are increased, 5r decreases and becomes negative 
while 5e increases. The vertical anisotropy is always negligi- 
ble with respect to the radial or azimuthal anisotropy (it is 
always negative in all performed simulations, as shown by 
the right panel in Fig. I19fl . 

Effect of resolution More resolved simulations are 
slightly more azimuthally anisotropic. As for some of the 
previously considered quantities, we observe that the com- 
ponents of the anisotropy vector, averaged in the vertical, 
azimuthal and in the usual radial range R £ [1,2.5], con- 
verge at A « 5-10^ (not shown), for the same reason already 
outlined several times. 

We note that radial anisotropy has been observed 
in MHD discs (see e .g. iFromang fc Papaloizoul 120061 : 
iFromang fc NelsonI l200q ) where the mechanism of merid- 
ional circulation is absent. Here we observe radially domi- 



nated discs only for low AV parameters, for which meridional 
circulation is also absent. 
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6.2 Power spectrum 

For each simulation we have computed the power spectrum 
of velocity fluctuations in a ring located at different radii 
from the central star and at several altitudes, as explained 
in Sect. [T2] As an example, in Figs. [20l and I2f I we present 
compensated spectra (divided by the Kolmogorov slope of 
—5/3, for easy comparison) evaluated at the midplane loca- 
tion [R, z)={2, 0) Different locations in the discs are consid- 
ered later on. 

Since we are dealing with an anisotropic system, we 
look at the one-dimensional power spectrum in the three 
directions: radial, azimuthal and vertical, which are given 
respectively in the first, second and third column of each 
row of the two figures. 

At least three regions can be identified in the power 
spectra extracted from simulations. These regions are de- 
fined by two wavenumbers: k\ < k2. The smallest scale 
corresponds to the resolution limit of the simulation that 
is given by the smoothing length at the considered position 
h{R, z) (since the discs are axisymmetric there is not depen- 
dence on 6) and defines the largest wavenumber ^2 (for each 
curve, this is marked by a circle in the plots). The region 
beyond fc2 corresponds to length scales inside the smooth- 
ing length, which are below the resolution. Therefore, the 
resolved region of the spectrum is to the left of the circle. 
The other scale marks the boundary between regions of the 
spectrum described by power laws P{k) — Pofc" with a dif- 
ferent slope: for k < ki the slope ai is close to zero, for 
ki < k < k2 the slope is called 02 and is always negative. 

In the case of isotropic turbulence one could identify 
fci with the 'forcing scale' and the region between fci and ^2 
with the 'inertial range'. However, in the present anisotropic 
case, it is not possible to find a direct link between fci and the 
forcing scale in the disc because of the aliasing phenomenon 
present in one-dimensional spectra of turbulent shear-stress 
fiuids (see e.g. |Popdl200Cl ). since the power corresponding 
to wavenumber fci also contains contribution from larger 
wavenumbers. 

In the following we focus on the qualitative evolution 
of the power spectrum with changing artificial viscosity and 
resolution. 

Effect of artificial viscosity. What is interesting is the 
dependence of the shape of the power spectrum on the AV 
parameters. Simulations with a < 1 and /3 < 2 show a pos- 
itive slope of the compensated spectra with a decrease just 
before the resolution scale {k < ^2), which mimics the decay 
corresponding to the dissipation scale. In this case a cascade 
is not clearly identified in any of the three directions. An ex- 
ample is given by simulation Sf with {a, j3) — (0.1,0) in the 
top row of Fig. [2S1 

For higher AV parameters the behaviour of the power 
spectrum depends on the component of the velocity fluctu- 
ations, confirming that we are dealing with an anisotropic 
system also at intermediate scales. In particular, for radial 
velocity fiuctuations in the case of q > 1 and /3 > 2 (e.g. 
middle and bottom left panels of Fig. I20|) and for the az- 
imuthal component when a ~ 5 and /? > 2 (bottom central 
panel of Fig. I20p , a progressively more extended fiat region 
(indicating a power spectrum with a slope close to the Kol- 
mogorov value) is visible between the positive slope region 
on the left and the decay immediately before the resolution 
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Figure 20. Compensated power spectrum k^i^P{k) of the veloc- 
ity field. Effect of changing a and /3. The three columns from left 
to right refer respectively to radial, azimuthal and vertical veloc- 
ity components. The power spectrum is computed along a ring 
located at _R = 2 and z = 0. The circle on each curve marks the 
wavenumber corresponding to the scale of the smoothing length. 
The number of particles is A^ = 2 ■ 10^ . 




Figure 21. Compensated power spectrum of the velocity field. 
Effect of changing changing the AV switch and the resolution. 
The circle on each curve marks the wavenumber corresponding 
to the scale of the smoothing length. 



scale on the right. The other azimuthal spectra show a pos- 
itive slope (e.g. top and middle central panels) , indicating a 
sub-Kolmogorov slope. Concerning the vertical component 
of the velocity, flat regions start to appear for values of a 
close to 1, just above the resolution limit. However, they are 
less extended than in the case of the other two components. 
For higher a, we observe that larger /3 leads to a steeper 
slope. 

In correspondence to the onset of meridional circulation, 
we observe the appearance of an approximately flat region 
in the compensated spectra for all velocity components and 
at the smallest resolved scales. This result suggests that the 



16 S. E. Arena and J.-F. Gonzalez 




Natural logarithm of the density 



Azimuthal acceleration 



Figure 22. Compensated power spectrum of the velocity field. 
Effect of changing the location inside the disc. Power spectrum 
computed along a ring located at i? = 2 and at three different 
heights above the midplanc for simulation S8 with (a, 13) = (1, 5). 



system is more isotropic at the small scales above the reso- 
lution limit. 

Effect of resolution In agreement with previous analy- 
ses, we observe a convergence of the power spectrum of ve- 
locity fluctuations when a number of particles A'' « 5 • lO'' is 
used. In addition, higher resolution discs are characterised 
by a cascade that is globally more extended in the wavenum- 
ber space than in lower resolution discs because they resolve 
smaller length scales, corresponding to larger wavenumbers. 
The resolution scale ^2 therefore shifts towards larger values 
(seeFig.[2lJ. 

In numerical simulations of turbulence a pile-up of en- 
ergy (positive slope of the spectrum) is often observed in 
the high wavenumber region of the power spectrum (cor- 
responding to the pre-dissipative region, close to grid res- 
olution) and is probably due to numerical dissipation. We 
observe a similar feature only in the power spectrum of ra- 
dial velocity fluctuations in intermediate and large viscosity 
discs (see left panels in Fig. [20l and [2T|) . The pile-up is re- 
duced by increasing resolution (see left panels in the top 
two rows in Fig. I2ip . supporting the numerical origin of the 
phenomenon. 

Simulations are observed to converge towards a power 
spectrum with a cascade close to the Kolmogorov one for 
the radial and azimuthal components. The system reflects 
its anisotropy in a cascade that is always more extended 
for radial fluctuations and very short for vertical velocity 
fluctuations. Such a particular feature of the cascade is due 
to the different physical scales present in the disc, which is 
more radially than vertically extended (the disc is thin). 

Changing location in the disc. The trends of the power 
spectrum observed with changing artificial viscosity and res- 
olution are qualitatively similar at different radial locations 
in the disc (not shown). In Fig. [22] spectra at different 
heights above the midplane are compared for each of the 
three velocity components in the case of simulation S8 where 
(a, /3) = (1, 5): moving away from the midplane, the slope of 
the cascade becomes slightly smaller and the cascade disap- 
pears due to decreasing resolution. In fact, at higher altitude 
the gas density is lower, implying a larger smoothing length 
(i.e. worse resolution). The wavenumber associated to the 
larger smoothing length is smaller and therefore it is closer 
to the wavenumber of the 'forcing scale', shortening the cas- 
cade. 

We conclude that the power spectrum preserves its 
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Figure 23. Probability distribution functions. Effect of AV pa- 
rameters. Left: natural logarithm of the density. Right: accelera- 
tion. PDFs are computed at the midplane and at distance R = 2 
from the central star. The density PDF is compared to a log- 
normal distribution (dashed line) in the top left panel. 



properties throughout the vertical and radial extension of 
the disc in the resolved regions of the wavenumber space. 



6.3 Intermittency 

In order to determine if intermittency is present in the sim- 
ulations we analysed the PDF of the density and of the 
acceleration in the azimuthal direction (the direction of the 
shear) and the 3^'^ and 4"' order moment of the PDF of 
several quantities. 



6.3.1 Density and azimuthal acceleration PDF 

Effect of artificial viscosity. The effects of the AV parame- 
ters on the shape of the PDF of the natural logarithm of the 
density are shown in the left panel of Fig. 1231 The distribu- 
tion of ln(p/(p)) is described by a Gaussian distribution. 
It follows that the density is characterised by a log-normal 
distribution, as observed in si mulations of supersonic com- 
pressible turbulence (see e.g. IPrice fc Federrathll2010l . and 
Sect. 1323]). 

However, in our simulations the standard deviation of 
the density logarithm ap, which controls the width of the 
distribution, is much smaller, since we are dealing with sub- 
sonic physical fluctuations (as shown in Sect. 15.31 Fig. [17] 
and Eq. I26|l . In addition, we observe that larger values of a 



Global flow and local fluctuations in 3D SPH discs 17 



Natural logarithm of the density Azimuthal acceleration 



Density 



Azimuthal acceleration 



1e+01 
1e+00 
1e-01 
1e-02 
1e+01 
1e+00 
1e-01 
1e-02 
1e+01 
1e+00 
1e-01 
1e-02 



(a,P)=(1,2) 



R=2 




(a,p)=(1,5) 




(a,P)=(1,0)^ 
Mu96 A 



N=5e4 
N=2e5 
N=1e6 





,.'-—•.,„ 



\, 



-0.2 0.2 

ln(p/<p>) 



-3-2-10123 

ae/<|ael> 



Figure 24. Probability distribution functions. Effect of resolu- 
tion. Left: natural logarithm of the density. Right: acceleration. 
PDFs are computed at the midplane and at distance R = 2 from 
the central star. The density PDF is compared to a log-normal 
distribution (dashed line) in the top left panel. 



and/or /3 lead to a narrower distribution. This behaviour is 
explained by the same reason: the Mach number of physical 
fluctuations decreases and leads to a smaller ap, which in 
turn produces a PDF close to a Gaussian distribution and 
gives values of the skewness and kurtosis close to Gaussian 
values (as we will show in Sect. 16.331 and as expected from 
Eq. I27p . The effect of changing /3 for a given a is smaller 
than that of changing a for a given /3, due to the fact that 
P controls the second order term in the IMG83| artificial vis- 
cosity. 

It should be noted here that even for high a and/or /3, 
although the effective viscosity becomes larger, the Reynolds 
stress is always non zero and still contributes to it. Therefore 
physical fluctuations are present and the disc is not laminar. 
For example, in our discs when a ~ 5 the effective viscosity 
increases to a2D ~ 0.1 and the Reynolds stress amounts 
to qrs ~ 10~^, with a corresponding turbulent Reynolds 
number ReT « 10*, see Sect. O and Figs. [12] and [H If the 
disc were laminar, one would not expect a distribution of 
densities but a single value, and therefore a delta function. 

In the right column of Fig. [23] we show the PDF of 
the acceleration in the azimuthal direction since it presents 
an interesting feature: the distribution shifts towards higher 
values when both a and /3 are increased. This is an effect of 
the onset of meridional circulation in discs with high enough 
artificial viscosity. In fact, the presence of meridional circu- 
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Figure 25. Higher order moments. Skewness (top) and kurtosis 
(bottom) of density distribution (left) and azimuthal acceleration 
distribution (right). 



lation imposes a net positive radial velocity (outward fiow) 
to the gas around the midplane that, combined with the 
disc rotation, implies a spiral-like fiow. The acceleration is 
thus characterised by a positive average azimuthal compo- 
nent (in contrast to a pure circular fiow where the acceler- 
ation is only radial). The peak location moves towards the 
value ag/{\ag\) — l for larger AV parameters, however it does 
not reach unity even for the larger AV combination con- 
sidered here (a, /3) = (5,10) since a fraction of the particles 
still presents a small but negative azimuthal component of 
the acceleration. Finally, the PDF of the radial and verti- 
cal components of the acceleration (not shown here) present 
the standard features of a pure circular fiow: both of them 
have a symmetric shape with the peak located around the 
non-zero Keplerian value or around zero, respectively. 

Effect of resolution . The effect of resolution is shown 
in Figure 1241 For all the combinations of AV parameters, 
more resolved discs are characterised by a narrower density 
PDF (left column), due to the reduction in the numerical 
noise. The azimuthal acceleration PDF (right column) is 
only slightly affected, with the peak shifted towards smaller 
acceleration for higher resolution, due to a correspondingly 
lower strength of meridional circulation (see Sect. [42]) . As 
already found in previous analysed quantities, a number of 
particles as large as TV ~ 5 ■ 10'' guarantees a good degree of 
convergence. 

In conclusion, the PDFs of density fiuctuations repro- 
duce those expected for non intermittent subsonic turbulent 
fiows, confirming that the discs are not laminar. 



6.3.2 Higher order moments: S and K 

Effect of artificial viscosity. For a given a, larger values of 
P result in values of the skewness of the density PDF closer 
to the Gaussian value 5* = (top left plot in Fig. I25|) . The 
same effect is observed when /3 is kept constant and a in- 
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creased. All other quantities show Gaussian values indepen- 
dently of the AV parameter used, an example is given by the 
azimuthal component of the acceleration (top right plot in 
Fig. [25]). Deviations of the kurtosis from the Gaussian value 
K = 2) axe qualitatively the s ame as for the skewness. 

The simulations with the lMuOa artificial viscosity show 
systematically l ower va lues of S and K with respect to simu- 
lations with the[MG83 artificial viscosity and wi th the same 
a and /3 = 0. The reas on is th at for the same a,|Mu9y discs 
are more viscous than IMGSSI discs since artificia l viscosity 
is applied to all particles in contrast to the lMG83| case . This 
trend is in agreement with that observed in IMG83I simu- 
lations with changing AV parameters: more viscous discs 
present lower S and K. 

Effect of resolution . The values of S and K only show 
small fluctuations when the resolution is increased (not 
shown), suggesting that convergence for these global quan- 
tities is present even for a number of particle as low as 
Ar«2- 10^ 

In conclusion, no significant deviation from a Gaussian 
distribution is observed: intermittency is not present. This 
result agrees with the intermediate turbulent Reynolds num- 
ber Rct proper to the present simulations (see Sect. 15.1) . In 
fact intermittency is usually ob served in very high Reynolds 
number flows (e.g. lFrisch|[l996l ) . 



7 DISCUSSION AND CONCLUSION 

We have presented the characterisation of the global flow 
and of the statistical properties of the fluctuations in the ve- 
locity and density fleld of the gas in SPH disc models, which 
are gaseous discs simulated by means of the SPH method. 

An essential tool for most of our results is the new 
method we have introduced for the determination of the 
effective viscosity of three dimensional axisymmetric discs. 
It consists in fitting the analytic expression for the ra dial ve - 
locity derived from two dimensional models (see e.g. lTL02h . 
depending both on R and z, to the vertical velocity profile 
extracted from the simulation at different radial positions in 
the disc. 

We have focused on the effects of changing the number 
of particles N and the values of the two AV parameters a 
and /3. The relevant results are summarised in the following 
points: 

(i) We have confirmed and quantified the numerical role 
of the number of particles A'^, which contributes to the nu- 
merical component of the SPH fluctuations (see its defini- 
tion in Sect.[T]). In fact, for all studied quantities, an increase 
in the resolution of the simulation (increase in the number 
of particles) has the standard numerical role of convergence 
towards the physical solution, and also leads to a more ex- 
tended cascade in the power spectrum of velocity fluctua- 
tions (smaller scales and therefore higher wavenumbers are 
resolved) . 

(ii) More significantly, we have found that the artificial 
viscosity (through the two AV parameters a and /3) con- 
tributes both to the numerical and to the physical compo- 
nent of the SPH fluctuations, in contrast to the behaviour 
oiN. 

There exists a relationship between artificial viscosity 



and the following three quantities: (a) physical fluctuations, 
quantifled by qrs, (b) the effective viscosity of the disc, 
quantifled by a2D and (c) the numerical noise. In the hy- 
pothesis of an ideal noise-free numerical scheme, we expect 
that at low [a, /?) the effective viscosity is due to turbulent 
viscosity, directly produced by physical fluctuations. At in- 
termediate (a, P) the effective viscosity has a contribution 
not only from physical fiuctuations but also from the AV 
term, which reproduces the effects of turbulence without di- 
rectly modeling eddies and vortices in the spirit of the ass- 
disc model. Finally, in the high (a, /3) range, physical fiuctu- 
ations are negligible and the effective viscosity is dominated 
by the AV term. However, numerical schemes are affected 
by numerical noise, which masks physical fluctuations. 

In our case, noise is related to the AV parameters: only 
for (a, /3) above a threshold is the noise reduced enough to 
allow physical fluctuations to emerge. We have identifled the 
threshold at (a, /3) ~ (1)2). In fact, in the discs simulated 
above this threshold, we have observed a sharp transition 
both of the global and of the local behaviour. 

• Concerning the global behaviour, we have found that 
with increasing artiflcial viscosity the global flow of the 
gaseous disc evolves from a chaotic radial velocity struc- 
ture to a more ordered meridional circulation pattern char- 
acterised by larger accretion rates and with the Reynolds 
number Reegf, associated to the estimated effective viscosity 
Q2D, below the turbulent limit for all simulations (however 
we observe that some of the high resolution simulations are 
approaching it). 

• Concerning the local behaviour, we observe that, in 
parallel to the onset of meridional circulation, the average 
Reynolds stress switches abruptly from positive to negative 
values and converges to a value that contributes at the 10 
per cent level to the effective viscosity (qrs ~ 0.1a2D), with 
an associated turbulent Reynolds number ReT around the 
turbulent limit. In particular, we have found that with in- 
creasing artiflcial viscosity, physical fluctuations have the 
following properties: 

- velocity fiuctuations become more and more subsonic 
and shift from radial to azimuthal anisotropy; 

- a cascade appears in the power spectrum of veloc- 
ity fluctuations, which tends towards a Kolmogorov-like 
spectrum, particularly for the radial component; 

- turbulent diffusion appears at the same time as the 
cascade and then the diffusion coefficient tends to de- 
crease. Diffusion tends to be damped by the onset of 
meridional circulation; 

- the presence of meridional circulation leads to a shift 
of the peak of the azimuthal acceleration PDF; 

- none of the studied SPH disc models presents fiuctu- 
ations characterised by intermittency. 

With these results we can give a first answer to the 
question raised in the Introduction: can SPH disc models 
correctly reproduce both the observed properties of PPD 
and the expected effect of turbulence? 

In the case of SPH disc models with the IMGSS artifi- 
cial viscosity we have found that for AV parameters above 
a threshold that approximately corresponds to the standard 
values (q ~ 1 and /3 ~ 2), the effective viscosity mainly rep- 
resents viscosity due to turbulence (of unknown origin) in 
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the spirit of the |SS73| Q-discs where intermediate and large 
scale turbulent eddies and vortices are not resolved. How- 
ever, at the smallest resolved scales, fluctuations with prop- 
erties similar to those present in turbulent flows have been 
identifled. A meridional circulation pattern and turbulent- 
like fluctuations coexist in SPH disc models. The former 
shapes the global scale of the flow, reproducing the observed 
mass accretion rate onto the star while the latter is present 
locally at the smallest resolved scales and contributes to only 
10 per cent of the effective viscosity of the disc. The resolved 
physical fluctuations reproduce the main effects expected in 
a turbulent 3D disc. In fact, the simulated discs are charac- 
terised by: a subsonic fluctuating velocity field presenting a 
turbulent-like diffusion and a small Kolmogorov-like cascade 
in the power spectrum of velocity fiuctuations at the small- 
est resolved scales (see e.g. simulations S7 and S8). These 
simulations present values of the Reynolds stress and of the 
diffusion coefficient that are very close to those observed in 
MHD accretion discs (where the source of turbulence is the 
MRI). For values of a and /3 below the threshold, the numer- 
ical noise dominates and masks physical fluctuations, these 
models should be avoided. 

A direct detection of turbulent eddies would require a 
reduction of the AV term in order to increase the contribu- 
tion of the Reynolds stress. However, we have shown that for 
the two AV implementations considered here the numerical 
noise dominates and impedes the growth of eddies in the 
velocity field when the levels of AV are below the thresh- 
old: in such cases the positive effects of the AV terms (par- 
ticularly of the /3 term) in correctly describing SPH par- 
ticles trajectories and in avoiding particle interpenetration 
are missing. Therefore, in order to reach higher Reynolds 
numbers it is necessary both to reduce the effective viscos- 
ity of the disc without increasing the numerical noise and to 
increase the resolution. This is not possible with the two arti- 
ficial viscosity implementation we are analysing in this work 
(since reducing the AV coefficients leads to an increase of 
the numerical noise). Recently developed artificial viscosity 
switc hes (e.g iMorris fc Monaghanl Il997l : ICuUen fc DehnenI 
[2OI0I) could help in this direction. However, it is still not 
possible to simulate fluids at the high Reynolds numbers 
(Recff 3> 10^) expected in astrophysical flows with the cur- 
rent tools and computer power. 

We should stress that simulating very high Reynolds 
number discs would be necessary if one wanted to investi- 
gate the possibility of the existence of pure hydrodynamics 
turbulence in 3D accretion discs. (We remind the reader that 
ID accretion discs are known to be stable with respect to 
hydrod ynamics instabU ities due to the Rayleigh criterion, 
see e.g. lArmitagdl2007l .) However, for our purpose of repro- 
ducing the effects of turbulence produced by an unknown 
source, the effective viscosity combined with a low and pos- 
sibly intermediate Reynols stress is enough. In fact, these 
models represent a starting point for the study of dust dy- 
namics, which is affected by physical fluctuations at small 
scales and meridional circulation at large scales. 

In conclusion, in SPH disc models, which a priori in- 
clude only the ingredients of star gravity and physical-Uke- 
viscosity, we have found that the artificial viscosity term, in 
addition to modeling a physical-like bulk and shear viscos- 
ity, can also play the role of an implicit turbulence model. 
In our SPH disc models, where we do not add any initial 



turbulent velocity field, the implicit turbulence model sus- 
tains and organises the random fluctuations present in the 
initial conditions of the disc. This result is in agreement 
with recent indications that a turbulence model is implicitly 
present in the SPH scheme used to simulate homogeneous 
and isot ropic turbulence with periodic boundary condition s 
(see e.g. lShi et all2012l : lMonaghanll201ll : lEllero et al]|2010l ). 
Here we have analysed the more complex case of systems 
without boundary conditions and where anisotropic turbu- 
lence is expected. 

The effects of the density and sound speed profiles on 
the global and statistical properties of the gas as well as 
those of the initial set up, of different SPH kernels and of dif- 
ferent artificial viscosity implementations will be addressed 
in more detail in a future work. 
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APPENDIX A: DIAGNOSTIC DETAILS 

Al The turbulent viscosity coefficient ut 

From the turbulent viscosity hypothesis (|Popell200d ). the 
Reynolds stresses are written in Cartesian coordinates as 
follows: 



(uiUj) = —2uTei 
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where en is the rate of strain: 
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and k the turbulent kinetic energy: 

k = l^i'^iUi). 

In classical accretion disc theory, only shear viscosity is rel- 
evant. Therefore the only non- vanishing co mponent of th e 
stress tensor is the R6 component (see e.g lLodatdl2008l ). 
which in cylindrical coordinates becomes: 
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Given the axisymmetry of the disc, the second term on the 
right hand side of Eg. IX4l vanishes and the turbulent viscos- 
ity coefficient becomes: 
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Tiie qrs coefficient 

The |SS73| coefficient corresponding to the viscosity of 
Sect. IXTl is defined by the relation vt — orsCs// that com- 
bined with Eq. IA5I gives: 
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the case of quasi- Keplenan discs we have: 
v-^ — -Rf2k, with Wk the Keplerian velocity and 
r2k = VGMr'^'^ the corresponding angular velocity. For 
thin discs, the approximation r ~ R holds, with r the radial 
spherical component and R the radial cylindrical compo- 
nent. 

Remembering that the scale height of the disc is related 
to the sound speed and to the Keplerian angular velocity by 
H — Cs/Hk, the Orb coefficient takes the form: 

(uRUe) r^k _ 2 {uruq) 
' RiQ.^)' ci ^3 cl 
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A3 The power spectrum, PDFs, 5* and K 

We consider a ring made of A'^g points and centred at the 
origin of the disc, with selected radius Ra and height Zs- 
The density and the smoothing length of each grid point are 
computed by means of an iterating procedure using the two 
coupled equations proper to the SPH scheme: 
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where the subscript a refers to the grid point and the sub- 
script b to its neighbours. 

Once the smoothing length is known, the value of the 
desired quantities (e.g. velocity components for the power 
spectrum, density and azimuthal acceleration for the PDFs 
and higher order moments) at each grid point are computed 
by the SPH smoothing technique. Now the one dimensional 
list of values is known and used to compute the power spec- 
trum of velocity components, the PDFs and the correspond- 
ing S and K coefficients. 

This procedure is applied to each simulation snapshot, 
then values are averaged in time (for 15.9 orbits at 100 au). 
For S and K the standard deviation of the distribution of 
values in time has also been considered. 



APPENDIX B: PROCEDURE FOR FITTING 
THE EFFECTIVE VISCOSITY 

The effective viscosity a2D in the disc is derived by fitting 
Eq. [15] to the data computed for a given snapshot of the 
simulation by means of the following three steps: 

(i) Computation of data from the selected simulation 
snapshot: radial profile E(_R) of the surface density, verti- 
cal profiles p{z) and vr{z) of the volumetric density and 
radial velocity at radial position R. 

(ii) Determination of the parameters present in Eq. 1151 
The surface density power law p at the selected location 
R is derived by fit of the surface density profile. The scale 
height H of the disc at R and _Ro is derived by fit of the 
relative vertical density profiles. Note that q is constant, 
since simulations are locally isothermal, _Ro is the length 
unit and Csq = Hq = H{Ro) since Ro — G = M = 1. 

(iii) Determination of a2D by fit to the vr data. 

All fits are performed by the general least square 
method, with x^ function defined by: 



x^-E 
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where {xi,yi) are the N data points from simulations, at 
the associated errors and f{xi,a) the function to be fitted, 
which depends on the parameter a. 



